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ABSTRACT 

We have studied dust evolution in a quiescent or turbulent protoplanetary 
disk by numerically solving coagulation equation for settling dust particles, us- 
ing the minimum mass solar nebular model. As a result, if we assume an ideally 
quiescent disk, the dust particles settle toward the disk midplane to form a grav- 
itationally unstable layer within 2 x 10 3 -4 x 10 4 yr at 1-30 AU, which is in good 
agreement with an analytic calculation by Nakagawa, Sekiya, & Hayashi (1986) 
although they did not take into account the particle size distribution explicitly. 
In an opposite extreme case of a globally turbulent disk, on the other hand, the 
dust particles fluctuate owing to turbulent motion of the gas and most particles 
become large enough to move inward very rapidly within 70-3 x 10 4 yr at 1-30 
AU, depending on the strength of turbulence. Our result suggests that global 
turbulent motion should cease for the planetesimal formation in protoplanetary 
disks. 

Subject headings: dust dynamics — planetary systems: formation — planetary 
systems: protoplanetary disks 



1. Introduction 

It is believed that particle settling and growth are important processes leading to the 
planet formation in protoplanetary disks. Observationally some evidences of dust size growth 
have been proposed based on dust continuum emission from the protoplanetary disks, such 
as smaller power-law indices of spectral energy distributions (SEDs) in sub-mm to cm wave- 
length bands, and fainter trapezium feature of 10/im silicate emission, compared with those 
of the interstellar dust grains (e.g., Beckwith & Sargent 1991; Miyake & Nakagawa 1993; 
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Beckwith et al. 2000; Kitamura et al. 2002; van Boekel et al. 2003; Przygodda et al. 2003; 
Wilner et al. 2005). Meanwhile, theoretically the initial stage of dust evolution when micron- 
sized interstellar dust particles grow into centimeter-sized particles with settling toward the 
disk midplane have been studied analytically and numerically; for example, Weidenschilling 
(1980) and Nakagawa et al. (1981, 1986)'s works in a quiescent disk, and Cuzzi et al. (1993) 
and Weidenschilling (1997, 2004) 's works in local turbulence induced by the shear between 
the dust layer and the gas near the disk midplane (see also Weidenschilling & Cuzzi 1993 
and references therein). 

In addition to the shear-induced turbulence, it is thought that there exists global tur- 
bulent motion in the protoplanetary disks, which is caused by thermal convective and/or 
magneto-rotational instabilities (e.g., Lin & Papaloizou 1980; Balbus & Hawley 1991). The 
turbulent motion of the gas is known to affect the dust evolution processes; turbulence in- 
duced relative motion increases the mutual collision rate, that is, the growth rate of the dust 
particles (e.g., Volk et al. 1980), turbulent mixing motion lets the particles move diffusively 
(e.g., Cuzzi et al. 1993), and turbulent eddies trap the dust particles (e.g., Klahr & Henning 
1997; Cuzzi et al. 2001; Johansen et al. 2004). On the other hand, the disk instabilities, 
namely, the existence of turbulent regions depend on the spatial and size distributions of 
the dust particles (e.g., Mizuno et al. 1988; Sano et al. 2000; Nomura 2004). Therefore, 
self-consistent treatment of the evolution of the dust particles and the turbulent regions (the 
disk instabilities) is needed in order to understand the very beginning of planet formation 
process in protoplanetary disks before the dust particles settle toward the disk midplane and 
form a dusty layer, which could lead to the planetesimal formation. 

Moreover, recent observations have been providing a huge amount of data of spectra and 
SEDs of dust continuum emission as well as molecular line emissions from protoplanetary 
disks. Theoretically reproducing these observational data have been also developed using 
detailed disk models (e.g., Kenyon & Hartmann 1987; Miyake & Nakagawa 1995; Chiang & 
Goldreich 1997; D'Alessio et al. 1998; Dullemond et al. 2001; van Zadelhoff et al. 2001; 
Aikawa et al. 2002; Nomura & Millar 2005). Although rather simple dust models have 
been used in those previous models, the dust properties affect the physical and chemical 
structure of the disks, and then the observable properties very much (e.g., D'Alessio et al. 
2001; Aikawa & Nomura 2005). Thus, we should carefully model the dust evolution in the 
disks and compare the model predictions with the observational data in order to interpret 
the observations and understand what is actually going on in protoplanetary disks. The 
SEDs of young stellar objects have been tried to be modeled by numerically simulating the 
dust size growth and settling processes in the disks (Suttner & Yorke 2001; Tanaka et al. 
2005; Dullemond & Dominik 2005). 
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In this paper, we have studied basic behavior of dust evolution in a quiescent or globally 
turbulent protoplanetary disk, especially focusing on time scales and growing dust size. It 
is done by numerically solving the coagulation equation, using a simple disk model, as a 
first step of understanding the initial stage of planet formation and modeling observational 
properties of protoplanetary disks. In the following section, we present the disk model as well 
as the basic equations for the vertical and radial motion and the coagulation of dust particles. 
In §3, we numerically calculate the dust size growth and settling toward the midplane in a 
quiescent disk, and compare the result with that obtained analytically by Nakagawa et al. 
(1986; hereafter NSH86). We also discuss the dust evolution in a globally turbulent disk in 
§4. Finally, the results are summarized in §5. 



2. Basic Equations and Models 
2.1. Disk Model 

As a disk model, we adopt the minimum mass solar nebular model (e.g., Safronov 1969; 
Hayashi 1981; Hayashi et al. 1985) in order to examine basic behavior of the dust evolution 
under a simple physical condition of the gas and compare the results with analytic calculation 
(see §3). In this model the gas surface density profile is given by 

S gas = 1.7 x 10 3 (i?/AU)- 3/2 g cm" 2 , (1) 

where R is the radial distance from the central star. The surface density of dust particles is 

where the solid density of a dust particle, p s , is set to be p s = 2 and 1 g cm -3 for R < 2.7 AU 
and R > 2.7 AU, respectively, taking into account the effect of water ice sublimation. We 
note that each of the mass and the surface density distribution is one of the most unknown 
factors in modeling protoplanetary disks. Numerical calculation of dust evolution using 
some different parameters for the disk mass and the power-low index of the surface density 
distribution as a function of the radial distance is performed by Tanaka et al. (2005). The 
dust and gas temperature is given by 

T = 280(i?/AU)" 1/2 K. (3) 

In this paper we treat the region where the vertical distance from the disk midplane, Z, is 
smaller than the disk scale height, 

H = (v^F/2)(c s /fi K ) = 0.0472(i?/AU) 5/4 AU, (4) 
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where c s = (8kT '/Vm^) 1 / 2 (m M = 2.34mH is the mean molecular mass and mn the mass 
of an atomic hydrogen) is the mean thermal velocity and f2x — (GM^/R 3 ) 1 ^ 2 (G is the 
gravitational constant and M* the mass of the central star) is the Keplerian frequency. And 
we simply assume that the profiles of the temperature (eq. [3]) and the gas density, 

Pgas = S gas /v^FiJ = 1.36 x 10- 9 (i?/AU)- 11/4 g cm" 3 , (5) 

are uniform in the vertical direction. The mass of the central star is assumed to be M* = 



2.2. Equations of Motion of Dust Particles and Gas 



The equations of motion of dust particles and gas in the disk are given by 

dU . , TT . GM S 
— = -Ap ga& (U-u)-- 

and 



-Ap gas (U-u)-—^R, (6) 



— = -A Pdnst {u -U)- —^R — , (7) 

az ti Pg as 

where U and u are the velocities of dust and gas particles in the inertial frame of reference, 
and pdust the spatial mass density of dust particles, p gas the gas pressure. The drag coefficient, 
A, is given by 

cjp s a for a < l g , 

3c s l g /2p s a 2 for a>/ 9 , 

following Epstein's and Stokes' law (Epstein 1924; Stokes 1851), respectively. The symbol 
a is the radius of the dust particles and l g is the mean free path of the gas particles, given 
by l g = m M /(<7 mo iPgas) = 1.44(i?/AU) n / 4 cm, where a mo \ = 2 x 1CT 15 cm 2 is the molecular 
cross section. The shape of the dust particles is simply assumed to be a compact sphere 
in this paper. We note that the dust shape (mass/area ratio) affects the dust evolution 
through the drag coefficient and the sticking rate (see next subsection), and many numerical 
studies have dealt with the effects of the dust shape, taking into account fractal structure of 
dust aggregates (e.g., Weidenschilling & Cuzzi 1993; Ossenkopf 1993; Weidenschilling 1997; 
Suttner & Yorke 2001; Dullemond & Dominik 2005). 

Now, we assume that the disk is axisymmetric and rotates around the central star 
at nearly Keplerian velocity, and set the dust and gas velocities relative to the Keplerian 
velocity, Vk(= RQk4>), as V = U — v K and v = u — v K - As far as the dust particles are small 
enough, the timescale for initial velocity of dust particles decaying due to the gas drag force 
is much shorter than the Keplerian time and the timescale of collision between dust particles 
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(e.g., Nakagawa et al. 1981). Thus, the mean motion of dust particles becomes steady soon 
in a quiescent disk (see also §4.1). By setting d/dt = in the equations of motion, we can 
derive the terminal velocities of the vertical and radial motion of dust particles as 

Vz-v z = -(S%/D g )Z, (9) 

where D g = Ap gas , and 

V R - v R = - rjv K and V R = — — — 2 2 rjv K , (10) 

L> + S l K Pgas + Pdust L> +M K 

where D = A(p g&s + p dust ) and 



^= -^^%r = m x 1(r3( w /2 



(see NSH86 for details). 



2.3. Coagulation Equation for Settling Particles 

We solve the following dispersed coagulation equation numerically for simulating the 
size growth of settling dust particles in the disk, according to Nakagawa et al. (1981) and 
Nakagawa & Kohno (1999); 

^ + jfiWzWm = -mMi) X>(*,J>(j) + Irn^Pii-jJMi-JMj), (12) 
where 

pm i+1 / 2 

ip(i) = / p(m)dm (13) 

is the mass density of the dust particles whose mass ranges from mj_i/ 2 to m i+ i/2, and rrii = 
(m,i-i/2 + mj+i/2)/2 for i — 1, • • • , n. The dust mass is binned into n intervals logarithmically 
as m i+ i/ 2 = £TOj_i/2, where rrii = (4n/3)p s af (the dust particles are assumed to have a 
shape of compact sphere), a\ = lpm, n = 320, and e = \/2 are adopted in this paper (see 
Appendix). Here, the total mass density of dust particles at a given position and time is 
given by 

rm n+1/2 ™ 
Pdust = / p(m)dm = 2_^<p(i). (14) 

Jm 1/2 i=1 



The second term of the left hand side of equation (12) shows the mass transport of dust 
particles in the vertical direction. Now, the mean vertical velocity of the gas is negligible 
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(y z ~ 0), so the dust particles settle toward the disk midplane. The mean settling velocity 
of dust particles with mass m 8 is derived as Vz(i) = —(p gas c s / psCi^Q^Z (if Oj < l g ) or 
— (3p gas cJ g /2p s a 2 )Q 2 i Z (if Oj > / 9 ) from equations (8) and (9). In a turbulent disk the dust 
particles are transported by turbulent mixing in addition (see §4.3). The mass transport of 
dust particles in the radial direction is not solved in this paper for simplicity. 

The symbol (3(i,j) is related to the sticking rate of two colliding dust particles, and 
given by 

f3(i,j) = n(ai + ajYSVpJmimj, (15) 

where we simply assume the sticking probability of p s — 1 in this paper. We note that 
the sticking probability will depend on size, relative velocities, chemical composition and/or 
shape of dust grains (e.g., Weidenschilling & Cuzzi 1993; Weidenschilling 2004; see also 
references therein). Lower probability will make the timescale of the dust evolution longer 
(e.g., Tanaka et al. 2005). Here, we neglect fragmentation of dust particles, which could 
occur if the particles become large and their relative velocities with small particles become 
high enough (e.g., Dullemond & Dominik 2005). As the relative velocity between the dust 
particles, SV, we take into account the thermal Brownian motion, 

5V B =( — + — 

V 7r / \rrii m,j/ 

where k is Boltzmann's constant, and the velocity differences in the vertical and radial 
directions, 8Vz = Vz(i)—Vz(j) and SVr = Vji(i)—Vji(j), which are derived from equations (9) 
and (10), respectively. The azimuthal velocity difference, SV^, has very weak size dependence 
as far as the dust particles are small (e.g., NSH86); hence we neglect SV^. We adopt the 
relative velocity of 5V = (5V£ + SV^ + SVr) 1 / 2 in a quiescent disk. In a turbulent disk 
we take into account the turbulence induced relative velocity, SV t , in addition as 5V = 
(5V£ + 5V% + 5V 2 + 5V t 2 fl 2 (see §4.3). 

In our numerical calculation the spatial grid in the vertical direction is taken equally 
spaced into 20 intervals within < Z < H . In addition, within the lowest interval of 
< Z < H/20, we take logarithmically spaced 17 sub-intervals as Z t+1 = 2Z t in order to 
resolve the region near the disk midplane in a quiescent disk. In a turbulent disk only equally 
spaced 20 intervals within < Z < H are used without sub-intervals because the dust size 
distributions are almost identical between the lowest vertical disk layers and the calculation 
of the diffusive mass transport due to the turbulent mixing is very time consuming if we use 
such small spatial sub-grids (see §4.3). 



(16) 
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3. Dust Evolution in a Quiescent Disk 

By numerically solving the coagulation equation for settling dust particles (eq. [12]), 
we obtain the dust size distributions at a given time and disk height at the Earth's (1AU), 
Jupiter's (5.2AU), or Neptune's (30AU) orbits. As an initial condition, we assume that the 
dust particles are well-mixed with the gas and have a radius of a certain value, a ini t. We 
adopt Oini t = 10, 20, and 60 /iin for R — 1, 5.2, and 30 AU, respectively, in order to compare 
the numerical results with the analytic calculation by NSH86. These values correspond to 
the wavelength of the peak emission at the local temperature and do not have particular 
physical meaning. The initial condition of the dust size distribution, however, does not affect 
the results very much. 

First, we compare our numerical result of the dust settling time with that by NSH86. 
Figure 1 shows the time evolution of spatial dust mass distribution obtained by our numerical 
calculation at the orbits of R = (a) 1AU, (b) 5.2AU, and (c) 30AU. The vertical axis 
represents the dust surface density from Z = to a characteristic height Z = (k = 
l,---,4) at R, S(Z < Zk) = f_z Pdust{R, Z)dZ, divided by the total dust surface density 
there, Edust (eq.[2]). The values, Zi, Z 2 , Z 3 , and Z 4 , at R =1, 5.2, and 30 AU are listed in 
Table 1. These characteristic heights are defined in NSH86 (see Fig.l of NSH86); to put it 
briefly, the vertical velocity, V z (eq. [9]), of dust particles dominates the radial velocity, V R 
(eq. [10]), above Z = Z 1 , and the gas density, p gas , is larger than the dust density, Pdust) 
above Z = Z 2 . The vertical velocity, Vz, dominates the radial velocity, Vr, again below 
Z = Z% where the dust density is high enough that the gas drag force hardly affects the 
radial motion of the dust particles. If most dust particles settle below Z = Z 4 , the dust 
layer becomes gravitationally unstable and could fragment into planetesimals. The dashed, 
dotted, dot-dashed, and solid lines in Figure 1 represent the dust surface density below 
Z = Zi, Z 2 , Z 3 , and Z4, respectively. As time goes on, the dust particles settle toward the 
disk midplane and more mass is included in the lower layer of the disk. In Figure 2 we plot 
the dust settling time at which 70% of the total dust mass settles below Z k (k — 1, • • • ,4) 
at the orbits of R = 1AU (squares), 5.2AU (circles), and 30AU (diamonds). Together with 
them, the dust settling time obtained by NSH86 is also plotted (triangles with solid lines). 
The figure shows that the numerical results are in good agreement with the analytic results, 
and most dust particles settle below Z = Z4, which leads to the formation of a gravitationally 
unstable dust layer, within about 2 x 10 3 , 6 x 10 3 , and 4 x 10 4 yrs at the orbits of R =1, 
5.2, and 30 AU, respectively, in a quiescent disk. 

Next, Figure 3 shows the resulting size distributions of mass density of dust particles, 
(p(i), normalized by pdust at (a) R = 1AU, t = 1 x 10 3 yr, (b) R = 1AU, t = 2 x 10 3 yr; (c) 
R = 5.2AU, t = 3 x 10 3 yr, (d) R = 5.2AU, t = 6 x 10 3 yr; (e) R = 30AU, t = 1 x 10 4 yr, 
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and (f) R — 30AU, t = 4 x 10 4 yr. The dashed, dotted, dot-dashed, and solid lines in each 
figure represent the size distributions at the characteristic heights, Z = Z 1 ,Z 2 ,Z 3 , and Z 4 , 
respectively. The time used in Figure 3 a, c, and e corresponds to t — t\ and that in Figure 36, 
d, and / corresponds to t — t±, where t k denotes the time at which 70% of the dust particles 
settle below the characteristic height, Z = Z k , at each orbit, R (see Fig. 2). Now, we can see 
from Figure 3 that the larger dust particles which have grown at the disk surface layer settle 
more rapidly toward the midplane with growing larger and larger. These processes lead to 
a bimodal size distribution near the midplane at the inner disk (e.g., Weidenschilling 1997). 
The gaps appear in Figure 3a and b around a — l g — 1 cm where the drag coefficient, A 
(eq.[8]), begins to follow Stokes' law, rather than Epstein's law. The drag coefficient always 
follows Epstein's law at R — 5.2 and 30AU in this model, where l g is much larger (l g — 1 x 10 2 
and 2 x 10 4 cm, respectively). The timescale of the dust size growth and settling is shorter 
at the inner disk where the effect of gravitational force of the central star is stronger. 

Finally, we compare our numerical result of the evolution of dust particle radius with 
that by NSH86. In Figure 4 we plot the resulting largest dust radii at Z — Z k (k — 1, • • • , 4) 
and R = 1AU (squares), 5.2AU (circles), and 30AU (diamonds). The largest dust radii, 
a max) at Z = Z k are obtained applying a criterion, i max = max{ i | <p(i) /pdust > 10~ 8 }, to 
the dust size distribution at t — t k (cf. thin solid lines in Fig. 3). The evolution of dust 
radii obtained by NSH86 is also plotted (triangles with solid lines). The figure shows that 
the numerical results are in good agreement with the analytic results within a factor of two, 
except at Z — Z\. The dust radii at Z — Z\ are larger in the numerical calculation because 
the relative velocity between dust particles due to the thermal Brownian motion, which 
works efficiently for small dust particles, is taken into account in the numerical calculation, 
but not in the analytic calculation in NSH86. Both numerical and analytic calculations show 
that the dust particles grow and their radii finally reach about 20, 7, and 1 cm at the orbits 
of R =1, 5.2, and 30 AU, respectively, just before they settle below Z = Z±. 

We note that the orbital decay is very little in a quiescent disk; during settling from 
Z = H to Z = Z A , the dust particles move radially by AR = 2.2 x 10~ 3 , 0.20, and 2.8AU at 
the orbits of R — 1, 5.2, and 30AU, respectively, according to NSH86. 

Here it should be commented that although we simply assume a totally quiescent disk in 
this section, in reality it is expected that the shear between the dust layer and the gas induces 
turbulence locally near the midplane as the dust particles settle (e.g., Cuzzi et al. 1993; Wei- 
denschilling & Cuzzi 1993; Weidenschilling 1997; Cuzzi & Weidenschilling 2005). Actually 
if we compute the Richardson number, J = —(dpd ust /dZ)(p gas + p^st^Q^Z (dV^/ dZ)~ 2 = 
Z(rjRp gas )~ 2 (p gas + Pdust) 3 (<9pdust/<9Z)~ 1 (e.g., Sekiya 1998; Chandrasekhar 1961), we can 
find that J < 0.25 and the shear instability will occur below Z = Z 2 at each orbit. In 
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such a turbulent layer the dust particles do not concentrate in the midplane and migrate 
inward very rapidly once they grow large enough as we will see in next section (although it 
is different in point that this is local turbulence). 

In conclusion, our numerical calculation have confirmed that NSH86's approximate 
treatment of the dust size growth and settling processes are appropriate for describing the 
dust settling time and the evolution of the largest dust size in an ideally quiescent disk, 
although they did not take into account the dust size distribution explicitly. This will be 
because most dust mass is included in the dust particles with the largest sizes, and the 
dust settling time is controlled by the largest dust particles which have the highest settling 
velocity. 



In this section we discuss the dust evolution in the disk in which global turbulent motion 
exists induced by, for example, thermal convective and/or magneto-rotational instabilities. 



First, we examine the vertical motion of one dust particle, not the mean motion which 
we have treated in the previous sections. The equation of motion of a dust particle (eq. [6]) 
in the vertical direction is written as 



The vertical velocity of the gas in a turbulent medium is generally given by u z = vT z + u' z , 
where the overline means a time average and the prime is a fluctuation due to the turbulent 
motion. Now we put the mean velocity to be W z = since we assume the hydrostatic 
equilibrium. As the component of turbulent fluctuation, we simply adopt an oscillating 
motion of u' z = vt exp(iutt) , which models the motion of the largest turbulent eddy with a 
velocity v t and a frequency uj t {uj t ~ v t /l t where l t is the eddy size). In this case, equation 
(17) has a general solution that consists of a mean motion part, Z(t), and a fluctuation part, 
Z'{t), caused by the turbulent motion of the gas, 



4. 



Dust Evolution in a Turbulent Disk 



4.1. 



Vertical Motion 




(17) 



Z(t) = Z(t) + Z'(t), 



(18) 



where 



Z{t) = d exp(Ait) + C 2 exp(A 2 t), 



(19) 
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_ D g v t exp[t(u t t-S)} _ , Dgut 

Z W - wo 2 , ,2^2 , n 2 ,211/2' ^- tan Q 2 2" ^ 

In equation (19) C\ and C 2 are integral constants, and 

Ai,2 = -±(D a ±y/D*-4Z%) 

f —-D fl , -^/A, for L> 9 > 2fi K , 
\ —(Dg/2) =p for D 9 <2fi K - 

If we assume that the eddy turn over frequency is equal to the Keplerian frequency, u t = Or, 
the fluctuation part (20) becomes Z'(t) « (v t /c s )H exp[i(£l K t — 7r/2)], which means that the 
turbulent gas motion forces the dust particle to continue to oscillate vertically with an 
amplitude (vt/c s )H and a frequency f^K, independent of the dust particle size (e.g., Landau 
et al. 1967). If we think a more realistic case, the turbulent velocity of the gas, u' z , will be 
modeled by a superposition of eddies with various sizes, velocities, and frequencies, which 
is often decomposed into Fourier components (e.g., Landau & Lifshitz 1959). In this case 
the term for fluctuating motion (20) is also given by a superposition of oscillations with 
various frequencies. Now, the mean motion part (19) shows the settling of the dust particle 
toward the disk midplane. From equation (18) we can derive the particle velocity, which also 
consists of a mean motion, V(t), and a fluctuation, V'(t), as 

V(t) = V(t) + V'(t), (23) 

where V(t) = dZ(t)/dt and V'(t) = dZ'(t)/dt. For a small particle which satisfies D g > 
2Q K , we obtain Z(t) ~ Z exp[— (Q^/D g )t\ (Z is the initial value of Z{t)) and V(t) ~ 

-(nl/D g )z(t). 

We note that in a quiescent disk in which u z = 0, a general solution of equation (17) is 
simply given by the mean motion part, Z(t), in equation (19). Therefore, the dust particles 
always settle toward the disk midplane; the motion of a large dust particle which satisfies 
D g < 2f2x is oscillation around Z — 0, damped (that is, the particle settles toward the 
disk midplane) with a timescale of 2/D g , while a smaller particle (D g > 2Qk) settles with a 
timescale of D g /Vt^ without oscillation (e.g., NSH86). 

4.2. Radial Motion 

Next, we will discuss the radial motion of the dust particles in a turbulent disk. When 
the dust particles are small enough, their motion is strongly coupled with the turbulent gas 



(21) 
(22) 
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motion. If the particles grow and reach a critical radius, they are released from the turbulent 
eddy trapping and migrate toward the central star (see e.g., Klahr & Henning 1997 for more 
detailed description of the dust particle motion in a turbulent eddy). The critical radius, 
a cr jt, is roughly estimated by comparing the friction time between the gas and dust particles, 
Tf = 1/Dg, with the turnover time of the largest turbulent eddy, r e dd y = and given by 

f c sPgas /p s Vt K for a < l 9 , 

mt I (3c s p gas / 9 /2p s fi K ) 1/2 for a>l g . 1 ) 

The critical radii at R =1, 5.2, and 30AU are a cr it = 32, 80, and 6 cm, respectively. When 
the dust radius reaches a cr it and the friction time, Tf, becomes as long as the eddy turnover 
time, r e( jdy, the radial velocity of the particle becomes the maximum, Vr ~ tjv-k = 5 x 10 3 cm 
s _1 (see eq. [10]), and the particle migrates inward very rapidly with the timescales of 
R/V R ~ 1 x 10 2 , 5 x 10 2 , and 3 x 10 3 yrs at R = 1, 5.2, and 30AU, respectively (e.g., Adachi 
et al. 1976; Weidenschilling 1977). 



4.3. Dust Size Growth in a Turbulent Disk 

Next, taking into account the properties of vertical and radial motion of dust particles 
mentioned in the previous subsections, we will numerically simulate the dust evolution in 
a turbulent disk by solving the coagulation equation (12). In this simulation we artificially 
remove the dust particles whose radii reach a crit as they migrate inward very rapidly. Nu- 
merical simulation including radial mass transport of the dust particles should be done in 
future (cf. Weidenschilling 2004). 

As we mentioned in §2.3 the relative velocity between the dust particles induced by 
microscopic motion of the turbulent gas, 5V t , is taken into account in this numerical cal- 
culation. As the relative velocity we adopt the approximate treatment by Weidenschilling 
(1984), which reproduces the result of Volk et al. (1980)'s analysis of the nonlinear response 
of a dust particles to the turbulent gas motion with a Kolmogorov spectrum. In addition, 
we use Mizuno et al. (1988)'s formula when the friction time is shorter than the turnover 
time of the smallest turbulent eddy (see also Markiewicz et al. 1991). The adopted relative 
velocity is 

6V t ={^^^ 1/2vt ^ T/i <r,<r fe0 , 

I \T fi + Tfj ) T kQ 111 T ks + T fi T kQ 111 T ks + T fj V * ' ^ S ' 

where Tf i is the friction time between the gas and dust particles with a radius a i: that is, 
Tf = l/Dg for a = di. The times of r feo (= r cd d y = 1/^k) an d T ks = Re" l ^ 2 T ko are the turnover 
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times of the largest (l t ) and smallest (Re 3 ^k) turbulent eddies, respectively. The Reynolds 
number, Re, is estimated as 

Re = vtk/v = ac s H/u = 2 x 10 U (R/ AU)- 3/2 a, (26) 

where we adopt the molecular viscosity of v — c s l g /3 = (c s m M )/(3<7 mo ip gas ) (e.g., Jeans 
1916). In this work we calculate the dust evolution in a weakly and strongly turbulent disk, 
in which we adopt a = 1CT 4 (v t = 10~ 2 c s and l t = 10~ 2 H) and a = 1CT 2 (v t = 10 _1 c s and 
l t = lO^if), respectively. 

As the mass transport of dust particles in the vertical direction, we take into account the 
transport due to turbulent mixing as we mentioned in §2.3. If we separate the particle mass 
density and the velocity into mean and fluctuating parts, the mass flux in the second term 
of the left hand side of equation (12) is described as Vz(i) ■ (p(i) = Vz(i) ■ (p(i) + V z (i) ■ f'(i), 
where the overline means a time average and the prime is a fluctuation due to the turbulent 
motion. As the mean vertical velocity we adopt Vz(i) = —(Q^/D 9 )Z (see §4.1). The second 
term of the right hand side of the equation, V z (i) • <p'(i), is the correlation of the fluctuations 
and treated as turbulent mixing, following the gradient diffusion hypothesis; V z (i) ■ <p'(i) = 
—D [dip(i)/dz], which works so as to diffusively uniform the mass density gradient of ip(i) 
(we omit the overline hereafter) in the vertical direction as the dust particles move around 
from eddy to eddy. Forthe diffusivity, we adopt D = Wt/j/(l+r//r c dd y ) = olc s H / '(1+JIk/ D g ) 
(e.g., Cuzzi et al. 1993; Weidenschilling 1997). The equation (12) is, therefore, solved by 
adopting 

Vz(i) ■ <p{i) = -^Z<p(i) - D ^ (27) 
for simulating the dust evolution in a turbulent disk. 

In Figure 5 we plot the resulting time evolution of the surface density of the dust particles 
whose radii reach a cr i t , S(a > a CT it) = f_ H dZ Yll=i ■ (*crit corresponds to a cri t), divided 
by the total dust surface density, Sdust (eq.[2]). As mentioned before, we have removed those 
large particles in the numerical simulation, taking into account the rapid radially inward 
migration. The figure shows that more than 70% of the total dust mass moves toward the 
central star very rapidly within about 70, 9 x 10 2 , and 1 x 10 4 yrs in a strongly turbulent disk 
(a = 10" 2 ; solid lines), while about 5 x 10 2 , 3 x 10 3 , and 3 x 10 4 yrs in a weakly turbulent 
disk (a = 10~ 4 ; dashed lines), at the orbits of R =1, 5.2, and 30 AU, respectively. The 
timescale of the dust size growth is shorter at the inner disk where the particle density is 
higher. 

Figure 6 shows the resulting size distributions of mass density of dust particles, tp(i), 
normalized by p d ust,o at (a) R = 1AU, t = 70 yr, (b) R = 5.2AU, t = 9 x 10 2 yr, and (c) 
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R = 30AU, t = 1 x 10 4 yr, in a strongly turbulent disk (a = 10 2 ). Figure 7 is the same as 
Figure 6, but in a weakly turbulent disk (a = 10~ 4 ), at (a) R = 1AU, t = 5 x 10 2 yr, (b) 
R = 5.2AU, t = 3 x 10 3 yr, and (c) R = 30AU, t = 3 x 10 4 yr. We note that the normalization 
factor in Fiugres 6 and 7 is different from that in Figure 3; Pdust,o used in Fiugres 6 and 7 
is the initial dust density, pdust.o = 4.2 x 10~ 3 p gas and 1.8 x 10~ 2 p gas for R < 2.7AU and 
R > 2.7AU, respectively, while pa U st used in Figure 3 is the dust density at a specific time 
and spatial position defined in equation (14). The thick solid, dashed, and dot-dashed lines 
in each figure represent the size distributions at Z = H, 0.5H, and 0.1H, respectively. 
The time used in the figures is when 70% of the dust particles grow large enough to migrate 
toward the central star very rapidly at each orbit, R (see Fig. 5). We can see from the figures 
that at each orbit, R, the size distributions of smaller dust particles are almost identical at 
each height, while those of larger particles are very different. In the strongly turbulent disk 
(a = 10~ 2 ; Fig. 6) the smaller particles have similar size distributions at each height mainly 
because the turbulence induced motion dominates the relative velocity and the dust size 
growth (e.g., Weidenschilling 1984) in almost all disk heights, Z; i.e., 5V ~ 5Vt independent 
of Z. Meanwhile, in the weakly turbulent disk (a = 10~ 2 ; Fig. 7) the differential vertical 
velocity, 5V Z , dominates the relative velocity at Z pa H where the gravitational force in the 
vertical direction is strong, while the turbulence induced relative velocity, 5V t , is dominant 
near the disk midplane. Therefore, the dust particles grow more rapidly at the disk surface, 
Z pa H, and small particles are replenished from lower disk layers via the turbulent mixing, 
which works so as to uniform tp(i). Consequently, the size distributions of smaller particles 
are not very different at each height also in the weakly turbulent disk. While the diffusive 
motion of turbulent mixing is strong enough to prevent the settling for the smaller particles, 
the larger particles cannot be sustained because of weak coupling with the gas and settle 
toward the disk midplane. So, the larger particles near the disk surface deplete as we can 
see from the figures. The depletion is more remarkable in the weakly turbulent disk (e.g., 
Dubrulle et al. 1995; Cuzzi et al. 1996; Cuzzi & Weidenschilling 2005). The gaps appear 
around a =(7, 4, and l.l)a: -1 / 2 pm at R — 1, 5.2, and 30AU, respectively, in Figures 6 and 7, 
owing to the discontinuities of the approximate treatment of the turbulence induced relative 
motion at r/ = Tk s (see eq. [25]). 

Our result that the most dust particles migrate toward the central star at a very short 
timescale suggests that global turbulent motion should cease for the planetesimal formation 
in protoplanetary disks. Unless the strength of turbulence is weak, the density of dust 
particles around the disk midplane will be low enough because of the turbulent stirring 
so that they cannot collisionally grow into planetesimals, as has been noted also by some 
previous works (e.g., Stepinski & Valageas 1996; Cuzzi & Zahnle 2004; Weidenschilling 2004; 
Cuzzi & Weidenschilling 2005). In a quiescent disk the dust particles will settle toward the 
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disk midplane and form a dust layer in a sufficiently short timescale as we have seen in the 
previous section. Afterwards the shear between the dust layer and the gas will cause a local 
turbulent motion near the midplane. Detailed analysis of the evolution of dust particles in 
such a turbulent dust layer, although it is beyond a scope of this work, is in progress for 
further understanding of the planetesimal formation process (e.g., Goldreich & Ward 1973; 
Cuzzi et al. 1993; Weidenschilling 1995; Sekiya 1998; Ishitsu & Sekiya 2003; Youdin k Shu 
2002; Weidenschilling 2003; Youdin & Chiang 2004). 

5. Summary 

We have investigated the dust size growth and settling toward the disk midplane in a 
quiescent or turbulent protoplanetary disk by numerically solving coagulation equation for 
settling dust particles. 

Our result shows that the dust particles settle toward the disk midplane to form a 
gravitationally unstable layer at a short timescale (2 x 10 3 -4 x 10 4 yr at R =1-30 AU) if 
we assume an ideally quiescent disk. The radii of the largest dust particles just before the 
formation of the unstable layer are 20-1 cm at 1-30 AU. The resulting settling time and 
evolution of the largest dust radius in our numerical simulation are in good agreement with 
those obtained by the analytic calculation in NSH86, although they did not take into account 
the dust size distribution explicitly. This is because most dust mass is included in the dust 
particles with the largest sizes, and these particles control the dust settling time. 

Also, we have discussed the dust evolution in an opposite extreme case of a globally 
turbulent disk to find that the dust particles are forced to fluctuate by turbulent motion of the 
gas, and grow to be large enough (32-6 cm at 1-30AU) to move inward very rapidly within 
a short timescale (70-3 x 10 4 yr at 1-30 AU). Thus, our result suggests that the disk should 
be quiescent or the global turbulent motion should cease before most mass of dust particles 
accrete onto the central star, in order to form planetesimals in protoplanetary disks. Self- 
consistent treatment of the evolution of the globally turbulent regions and the dust evolution 
processes is needed in future work. In addition, the dust evolution in a locally turbulent 
motion induced by the shear between the dust layer and the gas near the disk midplane 
should be investigated for further understanding of the planetesimal formation process. 

We would like to thank the referee, Dr. S.J. Weidenschilling for his comments which were 
greatly helpful in improving our paper. Also, we are grateful to Makoto Kohno for arranging 
the numerical code for coagulation of settling dust particles. This work is supported by 
"The 21st Century COE Program of Origin and Evolution of Planetary Systems" and the 
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A. Testing Numerical Solution of Coagulation Equation 

It is known that a numerical calculation of coagulation equation with inappropriate 
conditions, for example, coarse mass bins, causes some serious errors like an artificial accel- 
eration of coagulation that could lead to an artificial runaway (e.g., Ohtsuki et al. 1990; 
Wetherill 1990). Here we test our numerical solution of coagulation equation by compar- 
ing it with an analytic solution, using different numerical conditions. A linear kernel of 
mimj(3(i, j) = Qijrii + rrij) (Q is a constant) is used as a coalescence rate. In this case the 
coagulation equation with an initial condition of (/9(i)/(m i /m 1 ) 2 | t= o = Pduste~ m ^ mi has an 
analytic solution, 

_ ^stexp[- ?7 + K/m 1 )(2-e-_3] 7i[2K/mi)(1 _ 01/2] (M) 



K/mi) 2 (m l /m 1 )(l-e-") 1 / 2 

Pdustexpj-T? - K/mQtl - (1 - e~ r >) 1 / 2 } 2 } 
2vr 1 / 2 K/m 1 ) 3 / 2 (l-e^) 3 / 4 

where h(x) is the modified Bessel function and r\ = Qpdustt (Safronov 1963, 1969). 



(A2) 



In Figure 8 we compare the analytic solutions (A2) (solid lines) and the numerical 
solutions (crosses) which have an initial condition for the dust mass distributing at m\. The 
mass distributions divided by pdust at i] = 3, 6, 9, and 12 are plotted. Different numerical 
conditions are used in each figure; (a) e = a/2, (p m - m = 0, (b) e = \/2, (p m - m = 0, and 
(c) e = \/2, </? m i n = 10~ 21 pdust, where e is the intervals of the dust mass bins, TOj+i/2 = 
£mj_i/ 2 (see §2.3), and we prohibit the collisional coagulation (f3(i,j) = 0) if ip(i) < </? m i n 
or (p(j) < Win ( e -g- 5 Ohtsuki et al. 1990). In this paper we choose the conditions e = \/2 
and Win = 10~ 21 pdust) with which the numerical solution is unlikely to cause the artificial 
runaway. 
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Fig. 1. — The time evolution of spatial dust mass distribution obtained by our numerical 
calculation at the orbits of R — (a) 1AU, (b) 5.2AU, and (c) 30AU. The dashed, dotted, dot- 
dashed, and solid lines represent the dust surface density below the characteristic heights, 
Z = Zi, Z 2 , Z 3 , and Z 4 , respectively, divided by the total dust surface density, Sdust- The 
thin solid lines at E(Z < Zj)/S dust = 0.7 are used in order to estimate the dust settling time 
in Figure 2. 
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Fig. 2. — The dust settling time at the characteristic height Z = (i = l,---,4) and 
R = 1AU (squares), 5.2AU (circles), and 30AU (diamonds), obtained by our numerical 
calculation. The dust settling time by NSH86 is also plotted (triangles with solid lines). The 
numerical results are in good agreement with the analytic results, and most dust particles 
settle below Z = Z 4 , within about 2 x 10 3 , 5 x 10 3 , and 3 x 10 4 yrs at the orbits of R =1, 
5.2, and 30 AU, respectively, in a quiescent disk. 
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a [cm] a [cm] 

Fig. 3. — The size distributions of mass density of dust particles, tp(i), normalized by pdust 
at each orbit, R, and time, t, in a quiescent disk. The dashed, dotted, dot-dashed, and solid 
lines represent the size distributions at the characteristic heights Z = Z 1 ,Z 2 ,Z 3 , and Z4, 
respectively. The thin solid lines at yVPdust = 10~ 8 are used as a criterion to estimate the 
largest dust radii in Figure 4. 
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Fig. 4. — The largest dust radii at the characteristic height Z = Z^ (i = 1, •••,4) and 
R = 1AU (squares), 5.2AU (circles), and 30AU (diamonds), obtained by our numerical 
calculation. The evolution of dust radii by NSH86 is also plotted (triangles with solid lines). 
The numerical results are in good agreement with the analytic results within a factor of two, 
except at Z — Z 1 , and the dust radius finally reaches about 20, 7, and 1 cm at the orbits of 
R =1, 5.2, and 30 AU, respectively, just before most dust particles settle below Z = Z 4 . 
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Fig. 5. — The time evolution of the surface density of the dust particles which are larger 
than a crit in strongly (a = 1CT 2 ; solid lines) and weakly (a = 1CT 4 ; dashed lines) turbulent 
disks. Most mass of dust particles is included in large particles, which can migrate toward 
the central star, at a very short timescale (~ 70-3 x 10 4 yr at 1-30 AU). 
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Fig. 6. — The size distributions of mass density of dust particles, (p(i), normalized by Pdust,o 
at each orbit, R, and time, t in a strongly turbulent disk (a = 1CT 2 ). Note that the 
normalization factor is different from Figure 3 (see text). The thick solid, dashed, and dot- 
dashed lines represent the size distributions at Z = H, 0.5H, and 0.1H, respectively. The 
size distributions of smaller dust particles are almost identical at each height mainly because 
the turbulence induced motion dominates the relative velocity in almost all regions. Larger 
particles near the disk surface deplete since the turbulent diffusion against the settling toward 
the disk midplane is not strong due to weak coupling with the gas motion. 
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Fig. 7. — The same as Figure 6, but in a weakly turbulent disk (a = 10~ 4 ). The size 
distributions of smaller dust particles are not very different at each height because of the 
turbulent mixing in the vertical direction. The depletion of larger particles at the upper disk 
layer is more remarkable than that in the strongly turbulent disk. 
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Fig. 8. — Analytic (solid lines) and numerical (crosses) solutions for the mass distributions 
at T) — 3, 6, 9, and 12. Different numerical conditions are used in each figure; (a) e = \/2, 
V?min = 0, (b) e = \/2, </? min = 0, and (c) e = \^2, <p min = l(T 21 pdust- The conditions e = \/2 
and (fmm = 10~ 21 pdust are used in this paper. 
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